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Abstract 

An old conjecture in delay equations states that Wright's equation 

has a unique slowly oscillating periodic solution (SOPS) for every parameter 
value a > 7r/2. We reformulate this conjecture and we use a method called vali- 
dated continuation to rigorously compute a global continuous branch of SOPS of 
Wright's equation. Using this method, we show that a part of this branch does 
not have any fold point, partially answering the new reformulated conjecture. 

1 Introduction 

In 1955, Edward M. Wright considered the equation 

y'{t)^-ay{t-l)[\ + y{t)l a>Q, (1) 

because of its role in probability methods applied to the theory of distribution of prime 
numbers, and he proved the existence of bounded non constant solutions which do not 
tend to zero, for every a > 7r/2 [24]. Throughout this paper, we refer to equation (1) 
as Wright's equation. Since the work presented in [24], equation (1) has been studied 
by many mathematicians (e.g. see [4, 10, 11, 12, 13, 14, 19, 20, 21]). In 1962, G.S. 
Jones proved the existence of periodic solutions of (1) for a > Tr/2 [10]. Then in [11], 
he studied their quantitative properties and he made the following remark. 

The most important observable phenomenon resulting from these numeri- 
cal experiments is the apparently rapid convergence of solutions of (1) to 
a single cycle fixed periodic form which seems to be independent of the 
initial specification on [—1,0] to within translations. 

The cycle fixed periodic form he refers to is a slowly oscillating periodic solution. 

Definition 1.1. A slowly oscillating periodic solution {SOPS) of (1) is a periodic 
solution y{t) with the following property: there exist g > 1 and p> q+1 such that, 
up to a time translation, y{t) > on (0,(7), y{t) < on {q,p), and y{t + p) = y{t) for 
all t so that p is the minimal period of y{t). 
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After Jones made the above remark, the question of the uniqueness of SOPS in 

(I) became popular and is stiU under investigation after ahiiost fifty years. 

Conjecture 1.2. For every a> ^, (1) has a unique SOPS. 

It is worth mentioning that if Conjecture 1.2 is true, then the unique SOPS at- 
tracts a dense and open subset of the phase space (e.g. see [16]). Let us reformulate 
Conjecture 1.2, considering the partial work that was done since Jones's comment in 

[II] . In 1977, Chow and Mallet-Paret showed that there is a supercritical (forward 
in a) Hopf bifurcation of SOPS from the trivial solution ai a — ir/i [4]. We denote 
this branch of SOPS by J-q. In 1989, Regala proved a result that implies that there 
cannot be any secondary bifurcation from [22]. Hence, To is a regular curve in 
the {a,y) space. In 1991, Xie used asymptotic estimates for large a to prove that for 
a > 5.67, (1) has a unique SOPS up to a time translation [25, 26]. Here is a remark 
he made after he stated his result on p. 97 of his thesis [25]. 

The result here may be further sharpened. However, [. . .] the arguments 
here can not be used to prove the uniqueness result for SOPS of (1) when 
a is close to -j. 

Hence, his method might help to decrease the value 5.67, but new mathematical 
ideas are required to solve Conjecture 1.2. Based on the above discussion, here is a 
reformulation of the remaining parts of the conjecture. 

Conjecture 1.3. Denote by J-q the branch of SOPS that bifurcates (forward in a) at 
tt/2. Then 

1. J-Q does not have any fold in a ^ (^,5.67]; 

2. there are no connected components (isolas) of SOPS in a E (^,5.67]. 




Figure 1: Conjecture 1.3 fails if in the parameter range corresponding to a G 5.67], 
there exists a fold on or there exists an isola Ti of SOPS. 

In this paper, we propose to use a method called validated continuation in the 
parameter a to partially prove the first part of Conjecture 1.3. This method was 
originally introduced in [5] as a computationally efficient tool to compute equilibrium 
solutions of partial differential equations (PDEs) with polynomial nonlinearities. It 
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was then adapted to compute equilibria of PDEs for large (discrete) range of parameter 
values [7]. Afterward, it was combined with variational methods and tools from 
algebraic topology to prove the existence of chaos for a class of fourth order nonlinear 
ordinary differential equations [1]. In [2], validated continuation was generalized to 
compute global smooth branches of solution c;iirvcs of differential equations, both in 
the context of parameter and pseudo-arclength continuation. Finally, in a forthcoming 
work, the method is adjusted to compute equilibria of high dimensional PDEs [6] . In 
this paper, we use the theory of validated continuation developed in [2] to compute a 
global continuous curve of SOPS of Wright's equation. 

Theorem 1.4. Let e = 7.3165 x 10"''. Then the part of J^o corresponding to the 
parameter range a G [f + e, 2.3] does not have any fold. 



y 




Figure 2: Geometric representation of the result of Theorem 1.4. This curve represents 
a rigorous computation of a section of the set JFq. On the picture, the vertical axis is 
given by \\y\\ = sup{\y{t)\ ; t € [0,p], where p is the period of y}. 

For a geometric representation of Theorem 1.4, we refer to Figure 2. Before 
going into the details of the proof, let us make a few comments on the statement 
of Theorem 1.4. The reason why the result is valid only up to a = 2.3 does not have 
any theoretical justification. This is purely computational. In fact, when a grows, the 
proof becomes computationally difficult mainly because of the following facts. First 
of all, our computer-assited proof requires the computation of several sums which we 
compute using iterative loops with the Matlab interval arithmetic package Intlab [23] 
which is slow to evaluate loops of large size. A second observation is that the step 
size Act in the parameter a decreases significantly when one increases the parameter 
a. Hence, for larger a, the rigorous continuation still runs, but the step size decreases 
significantly. We come back to these issues in Section 6, where we make suggestions 
on how to possibly improve the result of Theorem 1.4. 

Another comment regarding Theorem 1.4 is that validated continuation in a can- 
not help ruling out the existence of a fold in the parameter range a G ]7r/2, 7r/2 + e[. 
This is due to the fact that the method requires having contractions which are uni- 
form in the parameter a. Because the trivial periodic solution y = is non hyperbolic 
at a = tt/2, the uniform contraction in the parameter a fails to exist near a = Tr/2. 
That raises the following question: How can we make sure that the global branch 
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of SOPS obtained with validated continuation for a e [7r/2 + e,2.3] actually comes 
from the Hopf bifurcation at a = 7r/2? It turns out that we can regularize the prob- 
lem at a = 7r/2 with the change of variable y{t) = l3z{t) and obtain a new problem 
(with continuation parameter (3 > 0) having a non trivial hyperbolic periodic solu- 
tion z{t) at 6 = and a = tt/2. This new problem, having now a as a variable 
(as opposed to a parameter), can be studied with validated continuation again, since 
uniform contractions can be proved to exist near (3 = and a = 7t/2. This is done in 
Section 5.4, where a rigorous continuation in the new parameter /3 > is performed 
in order to show that the branch of SOPS that we computed in the parameter interval 
a e [7r/2 + e, 2.3] is in fact the one that bifurcates from the trivial solution at a = 7r/2. 

Finally, it is important to mention that the value of e can be made smaller using 
our method. The choice of £ = 7.3165 x 10"* is made arbitrarily and we believe 
that with significant extra computational effort, this value can be pushed down up to 
e = 1 X 10~®. Once again, we discuss this possible improvement in Section 6. 

The paper is organized as follows. In Section 2, we transform the study of periodic 
solutions of (1) into the study of the solutions of a parameter dependent infinite 
dimensional problem f{x, a) = 0. In Section 3, the problem f{x, a) = is modified 
into an equivalent fixed point problem T{x,a) = x, whose fixed points correspond 
to zeros of /. The equivalence of the problem is shown and the functional analysis 
setting is introduced. In Section 4, we introduce the validated continuation method 
in the fashion of [2]. In Section 5, we prove Theorem 1.4 and finally, we conclude with 
possible improvements in Section 6. The computer programs used to assist the proof 
of Theorem 1.4 can be found at [9]. 



2 Set up of the problem f{x,a.) = 

The goal of this section is to transform the problem of looking for periodic solutions 
y{t + p) = y{t) of {1) into the study of the solutions of a parameter dependent infinite 
dimensional problem f{x,a) = 0. Let us introduce L to be the a priori unknown 
frequency of the periodic solution y. In other words, p = Hence, consider the 
following expansion of the periodic solution y in Fourier series 

oo 

y{t)= Cfce^'^", (2) 

k= — oo 

where the Ck arc complex numbers satisfying c-k = Ck- This is due to the fact that 
y G M. Plugging the two expressions 

oo oo 

y{t-l)= Cke-'^^e'^^* and y'{t)= ^ CkikLe'''^^ 

k= — oo k=—oo 

in (1) and putting all terms on one side of the equality, one gets a new problem to 
solve for, namely 

oo 
fco = — oo 

The left hand side of this last equation being a periodic solution with period one 
computes its Fourier coefficients by taking the inner product with e**^^*, for G Z. 



E 



ik\L ^ikiLt 



E 



ik2Lt 



= 0. 
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This procedure leads to the following countable system of equations 

gk = [ikL + ae-"'^] Ck + a ^ e-*'^^ ^Cfe, = 0, fc e Z. (3) 

Since c_fc = implies that g-k — Wki we only need to consider the cases fc > when 
solving for (3). Note that the frequency L oiy being unknown, we leave it variable and 
we are going to solve for it when solving / = 0. Denoting the real and the imaginary 
part of Ck respectively by Uk and bk, an equivalent expansion for (2) is given by 

oo 

y{t) = ao + 2 ^2 [c-k cos kLt — bk sin kLt] . (4) 

/c=l 

Note that ak = a-k and bk = —b-k- Hence, we get that bo = 0. Let 



Xk 



dof / (L, ao), k = 
{ak,bk), k>0 



and X = (xq, Xi, - ■ ■ ,Xk, ■ ■ ■ )"^. Let us denote by x^.i and x^. 2 the first and the second 
component of Xk € M^, respectively. In order to eliminate arbitrary shifts, we impose 
the normalizing condition y(0) = ao + 2 ^2^=1 <^fe = 0- Hence, let us introduce the 
following function h, which will ensure, by solving h = 0, that the scaling condition 
y{0) = is satisfied: 

00 

h{x) = ao + 2y^afc. 

fe=i 

For fc > 0, consider the real and the imaginary parts of gk, given respectively by 

Re{gk){x,a) = {acoskL)ak + {^kL + asinkL)bk (5) 
+a ^ {cos kiL){akiak2 - bk^bk^) + (sin /ciL)(afci6fc2 ^-bk^ak^), 

ki+k2 — k 

Im{gk){x,a) = — {—kL + a sin kL)ak + {a cos kL)bk (6) 
+a ^ - (sin kiL) {uk^ a^^ - bk^bk^) + (cos kiL) {ak^ bk^ +bkiak2)- 

ki+k2 = k 

Note that g-k = 'gk implies that Im{go) = 0. Hence, we do not incorporate Im{go) 
in the formulation of /. Hence, the function / is defined component-wise by 



fk{x,a) 



h{x) 
Re{go){x,a) 
Re{gk){x,a) 
Im{gk){x,a) 



, k = Q 
, fc>0 



Consider the notation fk,i (resp. fk^) to denote the first (resp. second) component 
of fk G K^. Defining / = {/fe}fe>o, we show in Section 3 that finding periodic 
solution y{f) of (1) satisfying y(0) = is equivalent to finding solutions of the infinite 
dimensional parameter dependent problem 

f{x,a)=Q. (7) 
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3 Set up of the fixed point equation T{x, a.) = x and 
functional analysis setting 

The purpose of this section is to transform the problem f{x,a) = into a fixed 
point equation T{x, a) = x. Then, the idea will be to apply an uniform contraction 
mapping argument on T. Let \is first put oursolf in a functional analysis setting by 
introducing a Banach space which is convenient for our study. The key ingredient in 
defining the space is that periodic solutions of Wright's equation are C°° [18]. This 
implies that the Fourier coefficients of the expansion (4) goes to zero faster than any 
algebraic decay. For s > 0, consider the weights 

= { \k\^ MO. 
These weights are used to define the norm 

\\x\\s = sup IXfelooWfe, (9) 
fc=0,l,... 

where \xk\oo = max{|xfe^i|, |a;fe^2|}, and the sequence space 

= {x = (a;o,a;i,a;2, . . .) , < oo}, 

consisting of sequences with algebraically decaying tails. Since the Fourier coefficients 
{xk}k>o decay faster than any given power of k, the set contains all sequences 
{L,ao,ai,bi, . . .) obtained from the Fourier expansion (4) of any periodic solutions of 
(1). We are ready to define the fixed point operator T. 

First of all, note that T will partially be constructed with the help of the computer. 
For that matter, we then truncate the infinite dimensional problem (7) into a finite 
dimensional one. More precisely, consider the finite dimensional projection /(™) : 
R^™ X R ^ M^™ defined component- wise by 

fk"'\xo,---,XTn.-i,a) = fki{xo,...,Xm-i,Ooc:),a) , k = 0,...,m-l, (10) 

where Ooo = (0)j>o- Consider a parameter value ao > tt/2. Recall from the discussion 
in Section 1 that since we aim for a contraction mapping argument, we consider only 
parameter values ao > tt/2. Indeed, at ao = '"'/S, the trivial solution is non hyper- 
bolic, meaning that Dxf{0,TT/2) is not injective. Suppose that at ao, we computed 
numerically x £ R^"* such that 

/M(5;,ao)«0. (11) 

This is done with a Newton-like iterative scheme. To simplify the presentation, we 
identify x = {L, ao, ai, bi,. . . , am-i, bm-i)^ with {x, Ooo)- Define 

A.-^(.,ao)-^g^^^^^^^ g^(,,ao)J- 

We use the subscript (•)p to denote the 2 (2m — 1) entries corresponding to A; = 
0, ■ • • , 2m— 2. Let Jp be a numerical approximation of the inverse of Da,/^^'""^) [x, ao), 
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02 be the 2x2 zero matrix and let 0^ be the 2 x 2(2rn — 1) zero matrix. Let 



J. 


0.^ 


0.^ 


0.^ 






02 


O2 


0. 


02 




O2 


0. 


02 


O2 


■'^2m+l 



(12) 



which acts as an approximate inverse of the Unear operator Dxf{x,ao)- More pre- 
cisely, given a; G O*, one has that 



Ax — i^JpXp , K2^_iX2m-li ■A2m'^2m) • • • ) • 

Lemma 3.1. Given (12) and (13), we have that A : $7* — *• 

Proof. First of all, there exists a constant 2x2 matrix S such that 



(13) 



|Afc-i| 



1„ 



for all k > 2m — 1 (see Lemma 5.3), where | • | means component- wise absolute values 
and <cw means component- wise inequalities. Considering a; e O*, one gets that 



\\Axl 



+1 



= max<^ |(yla;)o|oo, max \{Ax)k\^k''^'^ , sup \{Ax)k\ook''^^ 

y fe=l,...,2m-2 fe>2m-l 

= max<^ |(Jpa;p)o|oo, max \{JpXp)k\ook^^'^ , sup \h.k^ Xk\ook^~^^ \ 

I /c=l,...,2m-2 fe>2m-l J 

< maxi |(Jpa;^)o|oo,, max |(JpX^)fc|ooA;*+\ sup |Sa;fc|oofc' 

L fe=l,...,2m-2 fc>2m-l 



< 00, 



because ||a;||s = sup;^.>o |a;fe|ooWfe < 00 and S is a constant matrix. 



□ 



Let us comment on how, in practice, we make sure that the linear operator A is 
invertible. First of all, we verify that 



(14) 



with Ip being the 2(2rn — 1) x 2(2m — 1) identity matrix. If such inequality is satisfied, 

wc get that Jp is invertible. Recalling the definitions of fk^i and fk,i given in (5) and 
(6), respectively, and considering k > 2m — 1, we get that 



-Sk Tk 



(15) 



where Tk = aouo + ao(l + ao) cos kL and 6k = —kL + ao(l + ao) sin kL. Hence, a 
sufficient condition for to be invertible for all k > 2m — 1 is that 



m > 



ao|l + ao| 



+ 1 



(16) 
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Indeed, by (16), wo got that 6k < for all k > 2m — 1 and we can conclude that 
det{Ak) — + dl > 0, for all k > 2m - 1. Hence, if conditions (14) and (16) hold, 
the linear operator A defined in (12) is invertible. 

Given a parameter value a > ao, we define the fixed point operator T : x M to 
0« by 

T{x,a)=x-Af{x,a) (17) 

It is now important to remark that even if we constructed the operator T in a 
computer-assisted fashion, we still think of it as an abstract object. The finite part 
is stored on a computer, and the tail part, consisting of the sequence of matrices 
{A^^}fc>2m-i) is defined abstractly. 

Lemma 3.2. We have the following: 

(a) Let sq > 2 and fix a. Zeros of f{x,a), or, equivalently, fixed points ofT[x,a), 
that are in are in il* for all s > so- 

(b) Let s > 2. A sequence x ~ (xq, xi, a::2, . . . ) G is a zero of f, or a fixed point 
of T , if and only if y given by (4) is a periodic solution of (1) with y{0) = 0. 

Proof. For part (a), equivalence of zeros of / and fixed points of T is obvious, since 
the operator A is invertible. Suppose there exists x G 0*° such that f{x,a) = 0. 
Recalling that Xk = {ak,bk) for k > 1. that Ck — at + ibk and equation (3), we get 
that Qk = 0, for every fc > 0. Hence, for all fc > 0, we get that 



[ikL + 



ae 



ikL 



]ck = -a 



E 

fcl+fc2=fc 



(18) 



However, we have that 



J2 e-*'=^^Cfe,Cfe, 

fcl+fe2 = fc 



<2||x||^^ 



E 

fel+fc2 = fc 



< 



B 

k^' 



where B > is independent of k (see equation (38) in Lemma 5.2). Combining this 

inequality with (18), wc get that fc®"+'^Cfe is uniformly bounded. This implies that 
X G rj^o+i^ Repeating this argument, we can conclude that zeros of f{x,a) that are 
in are in for all s > sq. 

Finally, because the tail of a fixed point of T decays faster than any algebraic 
rate, all sums may be differentiated term by term, hence y defined by (4) is a periodic 
solution of (1) with y{0) = 0. On the other hand, any periodic solution of (1) is C°°, 
hence the tail of its Fourier transform decays faster than any algebraic rate, and thus, 
by standard arguments, the Fourier transform solves / = 0, and part (b) follows. □ 



We are now ready to introduce validated continuation. 



4 Validated Continuation 

Validated continuation [1, 2, 5, 6, 7] is a rigorous computational method to continue, as 
we move a parameter, the zeros of infinite dimensional parameter dependent problems. 
In our context, we use this technique to continue solutions of (7), as we move the 
parameter a. Lemma 3.2b shows that the problem of finding periodic solutions y of 
(1) such that 2/(0) = is equivalent to studying fixed points of T. We will find balls 
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in on which T, for fixed a, is a contraction mapping, thus leading to periodic 
solutions y of (1) satisfying y{0) = 0. 

Let Qo > 7''/2 considered in Section 3 and suppose that we computed a tangent 
X e R^" such that 



DJ^^\x,aQ)x+^—{x,ao) 
oa 



0. 



(19) 



As in Section 3, we identify x = {L, do, di, 6i, . . . , dm-i,i>m-i)^ with (x, Oqo)- Let us 
define the ball of radius r in (with norm || • , centered at the origin, 



B{r) n 



r r 



n 2 



so that a point b G i?(r) can be factored b = ur, with u G -B(l)- For 
we define the predictor based at ao by 



A„i; 



(20) 

a — ao > 0, 
(21) 



and balls centered at Xn 



B^^{r)=xa+B{r). (22) 

Definition 4.1. iei u,v € R"*^". VFe define the component-wise inequality by <cw 
and say that u <cw v if Uij < Vij, for all i = 1, . . . ,m and j = 1, . . . ,n. 

To show that T is a contraction mapping, we need component- wise positive bounds 
Yk = (pl'l^ , Zk = (^zl'l^ e for each k>0, such that, with = a - ao, 



[T{Xa, a) - Xa]k <cw ife(Aa), 



and 



sup 

6,ceS(r) 



[DxT{xa + b, a)c]k <cw Zk{r, Aa). 



(23) 



(24) 



Wc will find such bounds in Sections 5.1 and 5.2, respectively. Wc only consider 



A^ > 0, since we initiate the continuation at the parameter value ao 
move forward. The proof of the following Lemma can be found in [1]. 



e and 



Lemma 4.2. Fix s > 2 and a = ao + Aq,. If there exists an r > such that 

\\Y + Z\\s < r, with Y = (Yq, Yi, . . .) and Z = (Zq, Zi, . . .) the bounds as defined in 
(23) and (24), then there is a unique Xa G Bx^{r) such that f{xa,a) — 0. 

In order to verify the hypotheses of Lemma 4.2 in a computationally efficient 
way, we introduce the notion of radii polynomials. Namely, as will become clear in 
Sections 5.1 and 5.2, the fiuictions Yfe(Aa) and Zk{r,Aa) are polynomials in their 
independent variables. In fact, they are constructed to be monotone increasing in 
Aa- Also, for k > M = 2m — 1, where m is the dimension of the finite dimensional 
projection f^"^\ one may choose 



and 



Zk = Zm 



where ZM{r,Aa) >cw (o) is independent of k. The choice M = 2m — 1 will be 
justified in Section 5.1. This leads us to the following definition. 
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Definition 4.3. Let Yk{Aa) = (l) and Zk{r, Aa) = Znir, A^) for all k>M. 

We define the 2M radii polynomials {po, . . . ,pm-i,Pm} by 



Pk{r,Aa) = 



Yk{Aa) + Zu{r,Aa)-^^{\), fc = 0,...,M-l; 
ZM(r,A„)-^(;) k = M. 



The following result was first considered in [2]. 

Lemma 4.4. // there exists an r > and A^ > such that pk{r,Aa) < for all 
k = 0, . . . , M, then there exist a C°° function x : [ao, ao + A„] — > fl** : a x{a) 

such that f{x{a),a) = for all a G [ao,ao + A^]. Furthermore, these are the only 
solutions of f(x,a) = in the tube {a G [ao,ao + Aq,],^ — .t„ G B(r)}. 

Proof. By definition of the radii polynomials and because they satisfy Pk{r, Aa) < 
for all fc = 0, . . . , M, and by the choice of Yk and Zk for k > M, we get that 

\\Y + Z\U= sup ||n(A„) + Zfe(r,A„)||^u;fe <r. 

fe=0,l,... 

Since pk is increasing in A^ > (see Remark 5.5), existence and uniqueness of a 
solution x{a) for a G [ao, ao + A^] follows from Lemma 4.2. In particular, for every 
fixed a G [ao,ao + A„], T(-,a) : Bx„{r) Bx^{r) is a contraction. Consider the 
change of variable y = x — Xa- Then, the operator 

f : [ao, ao + A„] x B{r) ^ B{r) : (a, y) ^ f{a, y) = T{y + Xa,a) 

is a uniform contraction on B{r). Since / G C°° we have that T G 

C°° ([ao, ao + Aa] x B{r), B{r)). By the imiform contraction principle, we conclude 
that x{a) is a C°° function of a; sec e.g. [3]. □ 

The remaining part of the section is taken almost verbatim from [2] . 

In practice, we use an iterative procedure (with A^ varying) to find the approxi- 
mate maximal A^ (if it exists) for which there exists an r > such that the hypotheses 
of Lemma 4.4 are satisfied. If this step is successful, we let ai = ao + A^ and we ob- 
tained a continuum of zeros Co = {(a;''(a),a) | / (x°(a),a) =0, a G [ao,ai]}. We 
now want to repeat the argument with initial parameter value ai. Hence, we put our- 
self in the context of a continuation method, which involves a predictor and corrector 
step. Recalling the definition of the predictors based at ao given by (21), the predictor 
at the parameter value ai = ao + A^ is given by .ii = x + A'^x. The corrector step, 
based on a Newton-like iterative scheme on the projection /^™\ takes Xi as its input 
and produces, within a prescribed tolerance, a zero xi at ai. We can then compute 
a new tangent vector xi, built the new set of predictors xi + Aq.xi, construct the 
bounds Y, Z at the parameter value ai and try to verify the hypotheses of Lemma 4.4 
again. If we are successful in finding a new A° , we let a2 = ai -|- A^ and we get the 
existence of a continuum of zeros Ci = |(x^(a),a) | / (^x^{a).aj =0. a G [01,02]}. 
The question now is to determine whether or not Co and Ci connect at the parameter 
value ai to form a continuum of zeros Co U Ci . At the parameter value ai , we have 
two sets enclosing a unique zero namely 



Bo = xq + (ai - ao)i;o + -B(ro), 

and 

Si XI + Bin). 
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Figure 3: Bo H Bi contains a unique zero of (7) and Cq U Ci consists of a continuum 
of zeros. This picture illustrates the hypotheses of Proposition 4.5. 

We want to prove that the solutions in Bq and Bi are the same. We return now to the 
radii polynomials Pk{r, A^), k — 0, . . . , M constructed at basepoint (x, a) — {xi, ai), 
and evaluate them at Aq, — 0: 

Pk{r) =_pfc(r, 0). 

Since Pkifi) < 0, we find a non empty interval I = [r^, r^] containing ri such that 
Pk{r) are all strictly negative on I. We now have two additional sets enclosing a 
unique zero at parameter value ai , namely 

Bf x,+B{rt). 

The proof of the following result can be found also in [2] . 

Proposition 4.5. // Bq C B^ or B^ C Bq, then Co U Ci consists of a continuous 
branch of solutions of f{x,a) — 0, andCf)C\Ci ~ {(x"(ai), ai} — {(a;^(ai), ai} G 
BonSi. 

We have now all the ingredients to prove Theorem 1.4. 

5 The proof of Theorem 1.4 

The proof of Theorem 1.4 is constructive and it has two parts. The first one is a 
rigorous continuation in the parameter a G [7r/2 + £, 2.3] of a branch (denoted by Tq) 
of periodic solutions of (1). This part of the proof is presented in Section 5.3. The 
second part of the proof, presented in Section 5.4, verifies that Tq C J-q. In other 
words, we prove that the global solution curve .T-'q , computed in the first part, belongs 
to the branch of SOPS that bifurcates from the trivial solution at a = tt/2. 

Since we use validated continuation in the proof, we need to construct analytically 
the radii polynomials introduced in Definition 4.3. Section 5.1 is dedicated to the 
computation of the bound Y{Aa), defined component-wise by (23), while Section 5.2 
is dedicated to the computation of the bound Z{r,Aa), defined component-wise by 
(24). 
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5.1 The analytic bound Y{Aa) 

The goal of this section is to construct an analytic expression for the bound Y = 
Y{Aa) given by (23). Recall that this bound satisfies the following component- wise 
inequalities: 

[T{Xa,a) - Xa]k = [-Af{Xa.a)]^ <cwYk{Aa). 

As mentioned in Section 4, for a fixed value of ao, we consider a > ao and wc let 
Aq, = a—ao > 0. As a side remark, note that once the analytic bound Yk = Yk{Aa) = 
Yi.{q — ao) is dcrivc'cl. we use a computer program using interval arithmetic to get 
explicit numerical upper bound for 1^. By definition of fk given by (5) and (6), 
observe that fk{xa,a) = ( o ) for > 2to — 1. This is due to the fact that [xa]k ~ ( o ) 
for k > m. By definition of A given by (12), one can choose Yk{Aa) = (o), for 
k > 2m — 1. This fact justifies the choice of M = 2m — 1 already introduced in 
Section 4. Now that Yk is constructed for the cases k> M,we are left with the cases 
< < M — 1. Given i € {1, 2} and k € {0, . . . , 2m — 2}, let us compute the analytic 
bound Yh^i{Aa). As mentioned already in Section 4, we want to construct Yk^i{Aa) as 
a polynomial in Aq,. Recalling (23), we begin by splitting the expression f{xa,a) in 
two terms. The first term, very small because of the choices of x from (11) and x from 
(19), does not require any further analysis. The second term, not necessarily small, 
is expanded as an analytic polynomial using the software Maple and then bounded 
using further analysis. 

Let us now expand f{xa, cc) component- wise as powers of A^ using the function 

= fk,i{^c,a) = fk,i{x + {a- aQ)x,a). 

Recalling that A„ = a — ao > 0, Taylor's theorem implies the existence of ^ e 
[ao, a] such that 

fk,i{^a,a) = hl^iia) =hl^i{ao) + —^{ao){a-ao) + ^-T^{a*k,i)ia-aoy 



da 2 da^ 



d fk ' 

Dfk,i{x, ao)x + -^{x, ao) 



da 



Letting 



4°] = /fc,i(a;,ao); 
dg DfkAx, ao)x + %^{x, ao); (25) 

we have, as wanted, the following polynomial expression for fk^i, namely 

fk,i{x^, a) = dg + dg A, + dg(a^,,)A2 . (26) 

As mentioned above, the choice of the expansion (26) is made because the coefficients 
Ski and Ski from (25) are small. Indeed, d^.^] is small since {x,aQ) is a numerical 

approximation of (11) and <Sj}l is small because x is a numerical approximation of 

(19). In practice, d[". and d^,^] are evaluated using interval arithmetic. Hence, one can 
compute an explicit numerical upper bound for each of them. However, we cannot 
evaluate the quadratic coefficient d^^] (a^ ^) of (26) in the same fashion, because it 
depends on the unknown a^ ■ G [ao,a] = [ao, ao + A^]. The idea here is to define the 
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quantity aL'"'*^ = a^ . - ao e [0, A^] and to expand d)^/i{al i) = d^'^j-(ao + ^a'^') as 

powers of A^'^\ Once this expansion is done, the next step will be to use the fact 
that 

< A^f''^ < A„, for all i e {1, 2} and fc e {0, . . . , 2m, - 2}. (27) 

We will come back to (27) later. Using the mathematical software Maple, we compute 
analytic expressions d'j^l , dl^] , d[f] and d^j^l so that 



3(2), 



^k,i)^ 



dg(ao + A(^.^))=^dg (A(f.^)) 



i-2 



(28) 



j=2 



The Maple program D.mw generating the d):^, j = 2, 3, 4, 5 can be found at [9]. The 
first part of the program differentiate ^{a) = fk.i{xcnCt) twice with respect to a 

and then expands ~d^^\{ao + A^*^'*'') in powers of Ai'''*\ For more technical details 
about the expansion (28), we refer again to [9]. Combining (26) and (28), one gets 
that 

/.,(.„, a)=Y: dg Ai + dg (a2=.^))^" A^ 

As mentioned earlier, we now use property (27) and get rid of the dependence of 
fk,i{xa,oi) in terms of /sSa'^\ In order to do so, let us define 

= ((dg,d^]), (d«,dS), • • • , (d^£-2,i,d^^_,,,))'^ , i = 0, . . . ,5. 



For j = 2,3,4,5, let d^^^ 



) <i£f jO ) 



(a^^-)) 



J-2 



dC?) = ((dg^dgl),(dg,dg),...,(d^^_2,i,d2„_2,2 



and 



)) ,i = 



2,3,4,5. 



For the cases fc = 0, . . . , 2m — 2, we combine (27) and triangle inequality to obtain 
that 



\[T{Xa,a) -Xa\i 



-JFfF{Xa,a)\ 

j^J.d^K + EJ.d^P^l 

0=0 3=2 

1 5 



3=0 3=2 



Ai. 



As we mentioned before, the first part of the Maple program D.mw symbolically 
computes d^^-* , for j = 2, 3, 4, 5. The second part of D.mw helps obtaining the analytic 
upper bounds d[^'' (j 2,3,4,5) such that for i = 1,2, |d^-'|| < d[^\ The bounds 

£)g are presented in Table 1. It is important to note that all sums presented in 
Table 1 are finite sums. Hence, we can use a computer to compute them rigorously 
using interval arithmetic. Note also that Dq \ = for all j = 2, 3, 4, 5. Letting 



0,1 



|J^|DW , j = 2,3,4,5 
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k = 0, . . . , 2m — 2 




fclil (|i>fcl + Ibkl) + fciaj^ + ^^ok^L^b^ + ao'=i»fe - + »fc - ^"ofc^i-^a^ + kLb^ + <^ofciife 

+ E j-fcli ("fcj^ 3)12 ~ '''=1 ^k2) + ("fci ''*:2 + "ftl '''t2 + ''•'1 "'=2 + '"'^l "'=2) ^ i°o'»li^ ("fcj Sfcj + Sfcj Sj.^) 
/si +A;2 =fc 

-aofeli (Sfcj^ aj.^ + Sj.^ - Bj.^ Bj.^ - B^.^ Bj.^) + OQ (d^.^^ Bj.^ + ''fel "fe2)l 
+ E '=l^("fel^fe2 +^fel°fe2) + (»fel*fe2 +*fel»fe2 - ^fel''fe2 - ''fel^fe2) - i'»0feli^(afcia|c2 -^fel^fe2) 
+ <»0 ("fc, Afe, - Bfc, Bfc,) + aofeli (afc, Bfc^ + d^^ b^^ + b^^ d^^ + b^^ aj,,) 




2feidfc + ^aok^L^i,^ + ^k^LHk\ + 2fci,Bfc - ^aofe^i^dfc - ^fe^i^^^ 

_i_ \ 1 T . I 7i 1 -L/Si /ii _L 1 1 _L A 1 Pi 1 1 -1- ^f/Ti 1 -L T.lfit fit -1- A i /r i 1 

-1- 2^ ^ZKii^ yaf^^ Of^^ + »fc-L ''^2 + "fci "'fc2 ^ ^2 ^ ^ \ '=1 ^=2 ^1 ^2 ^ ^ -^OiO^l ^ ^"fc^ ^Jfeg + *'fci "fc2 } 

-\-k2 = k 

- ^aQkft^ (Sfci ^k2 + ^fc2 ^ ^^1 ^^2 ~ ^^1 ^^=2) ~ i'^l"^^ {^^1 ^^2 ~ ^^=1 ^^2) 1 

+ E -^''l^("fcl'^fc2 + '^?.-l^fc2 " ^fcl^fc2 ~ ^^l^/=2) + ^ ('^fcl^/^2 + ^A:i^fc2) " ^"0 ^ 1 ^ (^fc^ " ^i^fca) 
^1+^2='= 

-l«0'=?i'^(-,.T^>,., +-k, <^k;+H:, -A:, " ^^1^^^ k, ~^ k^, + ~^ k, - k-,,)\ 




^k^L^ (ibfcl + Idfcl) 

+ |3/=li (d^^ - ^ki^k2) + ■^'^O'^l-f'^ ^k2 + ^fcl "^^2) + h'^i^^ {^ki^k2 + ^ki^k2 + ^fci ^fc2 + ^fcl ^fc2)l 
fc ]^ + k2 = A; 

+ E \:ikiL {d^^b^^ + Bfcjdfc^) - ^aflfeli^ ("fel"fe2 " ''fel''fe2) ~ ^'"l^^ (°fel °fe2 + °fel °fe2 ~ ^fel ''fe2 ~ ''fel'fe2)l 




^fci^i^ [l-^fci^fca + ^fci "^^2 1 + l^^fci '^fc2 ~ ^fci ^^2 1] 



Table 1: The bounds D^' ■ 



we can finally set 



5.2 The analytic bound Z(r, A„) 

In this section, we construct analytically the bound Z = Z{r, Aq 
that this bound satisfies the component-wise inequalities 



(29) 



Recall from (24) 



sup [DxT{xa + b,a)c]k = sup [DxT{xa + ru,a)rv]k <cw Zk{r,Aa). 

b,ceB{r) u,veB{l) 



As mentioned previously in Section 4, we are going to construct each component 

Zk,i{r, Aa) (i = 1,2, fc > 0) of Z{r,Aa) as a polynomial in the variables r and 
Aa- In spirit, the construction of the polynomial expansion of Z(r, A^) is similar to 
the construction of the polynomial expansion of Y{Aa) of Section 5.1. We begin by 
splitting the expression DxT{xa + rw,, a)rv in two terms. The first term is small and 
does not require any further analysis. The second term, on the other hand, requires 
more analysis. It is expanded as an analytic polynomial using the software Maple and 
then bounded using analytic estimates. Let us now be more explicit. 
Introducing an almost inverse of the operator A defined in (12) 

I?x/('"-^nS,«0) 



0/^ 


0.^ 


0.^ 


^2m-l 

























we can split Df{xa + ru, a)rv into two pieces 

Dxf{xa + fu, Oi)rv = A^rv + [Dxf{xa + ru, a)rv — A^rv] . 
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Hence, we get 

D^T{xa + ru, a)rv = {[I - AA^]v) r-A [D^f{xa + ru, a) - A''] vr. (30) 

Note that the infinite dimensional vector [7— ^^''^jt; has only finitely many nonzero en- 
tries and its finite non trivial part, given by [/^ — J^D^f^'^"^^^'' (x, q:o)] e ]]j2(2m-i)_ 
has a small magnitude. This is due to the fact that is a numerical approximation of 
the inverse of Dj;f^'^"^~^\xp,ao)- In order to bound the second term of (30), further 
analysis is required. The idea is the following. First, expand each component of the 
term [D^fix a + ru, a) vr as a finite polynomial of the form 



([Z),/(s„ + ru, a) - At] vr) , . = ^ cg''^ V'^ A; 



Second, compute analytic upper bounds (7^'^''^^ so that \ci 



< C, 



(luh 



(uniform 



,(;i,;2)| 

with respect to i = 1,2). Finally, use the (7^'^''^^ to define the polynomial bound 
Z{r, A„). 

The computation of the c^'v'^^ is done analytically using the Maple program C.mw 
which can be found at [9] . The first part of the program computes an analytic repre- 
sentation of (Dxf{xa + ru,a)vr) I. ^. Then, ignoring the fact that the sin(-) and the 
cos(-) terms (coming from differentiating (5) and (6)) depend also on r and A^, it 
computes analytically, for all A; > and i G {1,2} the polynomial expansion 



([£'x/(a;„ + ru, a) - ^t] vr) ^ . = 



3 4-ii 



(31) 



Note that the coefficients c^'V'^'' i"^^) ^^iill depending on the sin(-) and the cos(-), 
which themselves depend on r and A^. The last part of C.mw is dedicated to the 
computation of the bounds Cj 



> such that 



--kA 



< Cj'"^^\ for i = 1,2. 



This part of the program uses several times the triangle inequality and the fact that 
cos I < 1. The bounds d^''^^^ are presented in Table 2. Note that the cases 



sm 



< C, 



(1,0) 
0,1 



IS 



C^'f^ and are treated differently. Indeed, the upper bound |cg^{"^| 

given in the first line of Table 2 and for the upper bound 1 4^^°^ I <C^'^2°\ we use the 
bound Cfe^'°^ (letting k = 0, this bounds is actually 0) on the second line of Table 2. 
Now that we have the bounds cj^^'''^\ we are ready to compute the bounds Zk{r, Aq,). 



5.2.1 The analytic bounds Zk(r, A„), fc e {0, . . . , M - 1} 



As mentioned earlier, the Maple program C.mw generates the coefhcients Cj. 
Defining C^J'M = : we get that 

/fe=0,...,M-l 

\[DxT{xa +ru,a)rv]^\ 

4 - J^D^f'^^'^-^\x,ao)] v^r-J^ [D^fixa+ru,a)rv-A''rv] , 



3 i-h 



3 i-h 



l = li2=0 
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fe = 0, t = 1 


(,(1.0) 


2 f ^ 
fe=2m — 1 '*^fc 


fc e {0, . . . , 2m - 2}. 


^(1.0) 


V- 1 ''-''1 1 1 '=-'«i 1 

fc^=2Tn — 1 '^'i^l 




oc kill kill 

fcl-oifc +fchodfc +Sfe| + ^ +°0'== |i| (hfcl + |5jc|) +2^5_li + + _L_i 

m-1 

+ H l^'l 1 l^fei^fe-ici + I'ki + "0 (»feiH-fei + "fei^fc-Jcj + 5jej dj,_;.j + Hi Sfe-fci)! 

fcl = -™ + l + fc 

771-1 

+ H l'=l 1 |-"fei »fe-l!i + 5fciEfe_)si - ao (<»fei Ofc-fei + -^ki''k-ki " 'fci^fe-lsi)! 

fcl = -77i + l + i: 

fcl=-m+l ^k-ki ki = -m + l + k 
^ m^l 2ao (Ifeil + |fc - fell) (|afc^ | + |6^^ |) 

ki = -m+l + k ^k-ki 


,(1.2, 


,n - 1 

'= (|^fe| + i^fci) + jz 1^1 1 l^fci ^fc-fci + "fci^fc-fci + ^fci "fc-fci + ^ki ^fc-fci + «o (^fci f'k-ki + ^fci '^fc-fci)! 

fcl = — m+l + fc 

fcj^ = -m+l+fc fcl = — m+1 '^k — ki 


,(1.3) 


fci = -m+l + fc 


,(2.0) 


^ ki+k2 = k^ki^k2 ki = -m+l '^k — k^ 

"^1 2<.o(|fcil + |fc-fcil)(|afei| + |EfcJ) , , 2n - - ,1 t 1,1-1 ,1 - h 

/ci = -77i + l + fc "fc-fci fci = -r7i+l + fc 


,(2.1) 


2.^ + E ^ + =(i^ii + i-^ii)fl''=^^"°'''=^l^l''=^^"''''=^l) 


,(2,2, 


e' 2(,.ii + ,.-.i„fKI + KI) 

ki=m~l \ '^k-ki / 


,(3,0) 


^ 4a;Q \ki 1 


,(3,1, 


E 



Table 2: The bounds c£'''^ for fc = 0, . . . , M - 1. 
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Before proceeding further, it is important to remark that the coefficients Cq^'°\ 
^(2,0)^^(2,1)^ ^(3,0) cj^'^^ of Table 2 involve infinite sums. This means that we 

have to use analytic estimates to bound these sums. The case of Cq''^^ is trivial. For 
instance, consider the estimate 



^ ±^ 1 



k=M 



(32) 



The infinite sums involved in cj^'^\cj^'^\ cj^'^^ and C'^'^^ can be bounded using 
the following result. 

Lemma 5.1. Let k G {0, . . . ,M — 1}, recall the definition of the weights Wk in (8) 
and define 



fe-i 



Then 



and 



E — 



fel+fc2=fc 



Ifcll - 1 



1 + 



< + — 4 + 



s-1 



1 + 



1^11 <- 

^ Wfc.Wfc, - {k + iy y ' S-2J ' 2'^" ' uJk ' (fc + l)«-i V ' s-1 

ki+k2=k 

Proof. First, 



k , k 
+ i;<Pk + — + - 



(33) 

(34) 

. (35) 



E — 



fe-1 



ki+k2=k 



y; 1 1 y-; 1 1^1 

4^ ojkiUJk-ki i^k kl{k - kiY UJk , tOki^Jk-ki 



< 



2 1 

^''^lJ^^^^^k{{k + ki) 



1 



— (4 



s-1 



Second, 



E 

ki+k2=k 



\ki\ 



^k-i_'^k2 



E 

A;i = — oo 



Ifcil 



^ki'^k-ki 



k-1 ^ 
E _ 



^ |fcl| 



1 



£ fcr'(fc + fci)' ' 2"'" ' oJk ' ^^{k + k^y-^ki 



k , k \ ^ 

o^^ + - + E 



< 



1 + 



+ TT'/'fc + — + 



(fc + 1)« V s - 2y 2^*" Wfe (A; + l)«-i V s-lj' 



1 + 



□ 



Hence, replacing the infinite sums of sums of Table 2 using the upper bounds (32), 

(34) and (35), we get new upper bounds C^'^''^-*. For k G {0, . . . , Af — 1}, we then 
define the Zk{r, A^) e to be the 2 dimensional fc*''— component of 

3 i-h 

Z,(r,A„) |[/, - J,i)/(2™-i)(5;,ao)] «.|r + ^ ^ | J, IC^'-'^^r'^ (36) 
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5.2.2 The analytic bound Zj^ir, Aq,) 



Consider k > M = 2m — 1. The goal of this section is to compute upper bounds 
CC'iM > such that for every k > M and i e {1, 2}, 



(37) 



where C^'^''^^ is independent of k and i. Wc computed the C^^'^) using the Maple 
program hatC.mw which can be found at [9] and by using the following result. 

Lemma 5.2. Defining 

s-2 



def rt 

7 = 2 



" M 


S 




+ 


_M -1 



41n(M - 2) _^ TT^ - 6 



M 



2 1 



and considering k > M, we have that 



Kl+fe2=K 



s- 1 



+ 7 



< 



and 



E 



ifci 



< 



1 



1 

M 



s-1 



+ 7 



fcl+K2=fe 

Proof. Let A; > M. By Lemma A. 2 in [1], we get 

fc-i 



2 7 
3 + 7 + 



s- 1 2 



(38) 
(39) 

(40) 



y t 



~ k^ \ 



< — 7- 



" k ' 




1. 


+ 



41n(A:-2) _^ tt^ - 6 



2 1 
k^2 



The rest of the proof is a minor modification of the proof of Lemma 5.1. 



□ 



The bounds (39) and (40) are used to find the C^'^''^) satisfying (37). The bounds 
are presented in Table 3. We still need one last estimate before defining the 
bound ZM{r, A^). 

Lemma 5.3. Let Z > 0, ao € K and consider m such that (16) is satisfied. Define 

M 



and 



ML-ao\l+do\ 
77-1 l7"o(|ao| + |l + ao|) 



> 



\ P M 

Then for all k > M, Ak is invertible and 

1, 



^ao {\ao\ + |1 + ciol 



k 



(41) 
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c(i.o) 



c(i.i) 



fci=i ' 



^7^ — ^ 



^ + (2ao|l + «ol + 1)1^1 + ^'1 + °^'^' +"1: -^(la^, +aodi, l + |Efe +aoi'feil) 1 



fcl=l ^ 



(i- Jji V 



C(1.3) 



C(2,0) 



+ E 2ao(|afcj| + |5fc^|) 
fcl=l 



2(1 + |do + aodol) + (-4 + 7^ + t) + E /"°' (|5fc^ + Q0<ifci I + lEfci + aof-fci I) I 1 + 



C(2.1) 



(l ~ 5m- 1) 



+ Yl 2(|aj,^ + aodfcj^ I + Ib^^ + agik^ I) 
fci =1 



2|i.oi + E :^ — ^<i"fcii + I'*.!!) 



7 I + 13 2(|afc | + |5fc I) 



1,1 ~ 2m-l j ; 



3(3,0) 



4oo (3 + ^ + ^) 



C(3.1) 



12 + + 27 



Table 3: The bounds (7(^1 



Proof. The fact that Afe given by (15) is invertible for all k > M > m follows from 
the choice of m given by (16) and we then get that 



1 



+ Sl 



Sk Tk 



Since k> M> 



\Sk 



Therefore, 



and then 



= kL — ao{l + ao) sin kL 
> fcL — ao 1 1 + flo I 

«o|l + aoT 
k 

ao|l + QqI 
M 

P 



= k\^L 

> k( L 



- > 0. 

P 



\Sk\ ^ k 



5k 



rl + SI 



14] ^ 

-'fe 14 1 



Finally, since |rfc| < ao (|ao| + |1 + ao|), wc get that 

Tk 



^ gp (|ao| + |1 + oqI) ^ gp (|ao| + |1 + ao|) 



< 



^1 + ^^ 
p2ap(|ap| + |l + ap|) ^ 1 



A;2 



A; 



SI 

p^ap(|ap| + |l + ap|) 
M 
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□ 



Wc arc now ready to define ZM{r, Aq) in the fashion of Definition 4.3. 
Lemma 5.4. Define 

o2 



and consider k> M. Then 



ao {\ao\ + |1 + ao\)+p 



3 4-h 

ii=i;2=o 



(42) 



/ M 

\[DT{x„+ru,a)rv]^^\ <cn, ZM{r,A„) ( — 
Proof. Let k > M. Combining equations (30) and (31), and Lemma 5.3, we get that 



[DT{xa + ru,a)rv]^ 



— ^ [^Df{xa + ru, a)rv — A^rv] ^ | 

3 4-(i 



< 



1 = 1/2=0 

3 4-;i 



< 



cw 



il = U2=0 



ZM{r,A^)i^—j . 



□ 



Remark 5.5. Recalling the definiUons ofYp, and Zm, given respectively by (29), 
(36) and (42), one easily observe that the radii polynomials pk(r, Aa) from Defini- 
tion 4-3 are monotone increasing in the variable A^ > 0. 



5.3 First part of the proof of Theorem 1.4: Rigorous compu- 
tation of the branch using vahdated continuation 

In Sections 5.1 and 5.2. we constructed the bounds Y and Z, respectively. The 
coefficients in Tables 1, 2 and 3 provide us an analytical representation of the radii 
polynomials associated to (7). The following Procedure is an algorithm to compute a 
global continuous branch of solutions of (7). 

Procedure 5.6. To check the hypotheses of Lemma 4-4 df^d Propostion 4-5 on the 
interval a £ [n/2 + s, 2.3], we proceed as follows. 

1. Consider minimum and maximum step-sizes Amin = 1 x 10~^^ and Amax = 2, 

respectively. Initiate s = 3, m = 6, M = 2m — 1, ao = tt/2 + e, ro = 0, 
Aq, = 5 X 10~^ G [Amin, Amax], A° = 0, and an approximate solution Xp of 
f^"^\-,ao) = given in Figure 4- Initiate Bq = Bx^{ro). 

2. With a classical Newton iteration, find near Xp an approximate solution Xp of 
f^"^\xp, ao) = 0. Calculate an approximate solution Xp of Dxf'^™'\xp, ao)Xp + 
Daf^"^\xp,ao) = 0. Using interval arithmetic, verify that conditions (14) and 
(16) are satisfied (this guarantees that the linear operator A defined in (12) is 
invertible). 
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3. Compute, using interval arithmetic, the coefficients of the radii polynomials pk, 
k = 0^ . . . , M given in Definition 4- 3. This is the computationally most expensive 
step, since it involves computing all coefficients in Tables 1, 2 and 3, and in 
particular requires the calculation of many loop terms. 

4- Calculate numerically! = [ri,ri] = PlftLoi^ — 0|Pfe(^)0) ^ 0}- Consider 
S" = Bs^ (rf) and = B^^ {rt). Verify that Bq c or B^ c Bq. 

5. Calculate numerically I = [/_,/+] = flfeloi'' — 0|Pfe(^;^a) ^ 0}- 

• If I = then go to Step 7. 

• If I ^ then let r = -j^- Compute with interval arithmetic Pk{r, Aa). 
If Pk{r, Aa) < for all k = 0, . . . ,M then go to Step 6; else go to Step 7. 

6. Update A° <— and r^ ^ r. If ^A^ < Amax then update A„ <— ^Aq, and 

go to Step 5; else go to Step 8. 

7. If A*^ > then go to Step 8; else if jqA^ > Amin then update A^ <— ^A^ and 
go to Step 5; else go to Step 9. 

8. The continuation step has succeeded. Store, for future reference, Xp, x^, tq, 
ao and A|^. Determine ai approximately equal to, but interval arithmetically 
less than, ag + A" . Make the updates ao <— ai, Aa ^ A°, x^, ^ x^ + A^i^ 
and A° ^0. // one of the last two components of Xp has magnitude larger 
than 1 X 10^^, update x^ ^ {xp,0,0), m ^ m + 1 and M 2m — 1. Update 
Bq ^ Bxp {ro) and go to Step 2 for the next continuation step. 

9. The continuation step has failed. Either decrease Amin and return to Step 7; or 
increase M and return to Step 3; or increase m and return to Step 2. Alterna- 
tively, terminate the procedure unsuccessfully at a = ao (although with success 
on [7r/2 -\-e,ao]). 



L 


1.570599180042083 


ao 





ai 


0.000393777377493 


bi 


0.031377227341359 


0,2 


-0.000389051487791 


b2 


0.000206800585095 


0.3 


-0.000004694294098 


b3 


-0.000001372932742 


04 


-0.0000000311811.38 


64 


-0.00()0()0()3."')0."')2(i66 


0-5 


-0.000000000114467 


65 


-0.000000000397361 



Figure 4: Approximate zero Xp at the parameter value ao = f + £■ 

The Matlab program intvalWrightCont.m, which can be foimd at [9], performs 
Procedure 5.6 successfully on the parameter interval [7r/2 + e,2.3]. Hence, by con- 
struction, we get the existence of a continuous one dimensional branch of periodic 
solutions To which docs not have any fold in the range of parameter [7r/2 + £,2.3]. 
This result follows from the uniform contraction principle and Proposition 4.5. The 
last step of the proof is to show that is the branch of SOPS of Wright's equation 
that bifurcates from the trivial solution at a = 1^/2. 
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5.4 Second pcirt of the proof of Theorem 1.4: Bifurcation anal- 
ysis at a = 7r/2 to show that ^* C 

In this section, we show that the branch J^q comes from the Hopf bifurcation at 

a = tt/2. For a detailed analysis of this Hopf bifurcation, wc refer to Section 11.4 of 
[8]. Consider the change of variable y{t) = Pz{t). Plugging y{t) = (3z{t) in Wright's 
equation (1), we get 

z(t) = -az(t - + (3z{t)]. (43) 

Consider the problem of looking for periodic solutions of (43), with the parameter now 
being /? > (a is now considered as a variable). We impose to the periodic solutions 
the conditions ^;(0) = and i(0) = — 1. More precisely, we consider the problem 



z{t) = -az{t - + 0z{t)], (3>0, 
z{t+^)=z{t), 
z(0) = 0, i(0) = -l. 



z(t+^)=z(t), (44) 



When /? = 0, q; = tt/2 and L = tt/2, equation (44) has solution z{t) = —-^ sin (f i). 
This solution corresponds to the Hopf bifurcation point y{t) = (— ^ sin (ft)) = 0, 
when a = tt/2 and L = tt/2. The idea is to use validated continuation (in the 
parameter /? > 0) on problem (44) and to connect the rigorously computed branch of 
SOPS of (44) to the left point of J^q. It is important to note that this new validated 
continuation cannot help ruling out the existence of fold in the space {a,y), but only 
in the space (/3, z). 

Considering the periodic solution z{t) in Fourier expansion, we do as in Section 2 
and consider a function to solve for. Defining X = {a,x), it can be shown that an 
equivalent problem of (44) is F{X,(3) = 0, where 



where 



and 



-l + 2LJ2Zikbk, k = -l 
a{ao + 0al + 2pj:Z=i(^osk,L){al+bl)) ' " ^ 



Rk{L,a) = 



k\ +k2 = k 



a COS kL —kL + asinkL 
a sin kL a cos kL 



(45) 



cosfciL sinfciL 
— sinfciL coskiL 



To apply validated continuation on problem (45), with (3 > being the parameter, 
wc need to construct the radii polynomials. Here, wc do not provide analytically the 
coefficients of the radii polynomials associated to (45), since they are similar to the 
ones associated to (7). A procedure similar to Procedure 5.6 is applied on (45) to 
get the existence of a continuous branch of SOPS of (44) on the parameter range 
P e [0,/3o], where Po = 0.099847913753516. We denote this branch by G^. See 
Figure 5 for a geometric representation of Go- the right most point of Gq, we have 
a set Bq containing a unique solution of F{X,(3a) = 0. Using a similar argument 
than the one presented in Proposition 4.5, we can show, via the change of coordinates 
y = (3z, that the solution in the set Bq and the solution on the left most part of the 
branch are the same. Hence, we proved that J^q G Tq. □ 
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'^0 GX12 004 005 OQS 0.1 



Figure 5: A branch of SOPS of (44) on [0,/3o]. 

6 Future Work and Acknowledgments 

As mentioned in Section 1, we believe that Theorem 1.4 could be improved signif- 
icantly. The reason why the proof was stopped at a — 2.3 is due to the fact that 
the Matlab program intvalWrightCont.m [9] becomes slow for large a. Indeed, the 
evaluation of the coefRcients of the radii polynomials is computationally expensive, 
mainly because of all the iterative loop evaluations in Step 3 of Procedure 5.6, a task 
that the interval arithmetic Intval is not efficient at doing. Using a different program- 
ming language (like C or C + +) would decrease significantly the computational time. 
We believe that we could push the parameter value up to a = 3 using a C program. 
This speculation is based on simulations that were done in Matlab without interval 
arithmetic. We could, with the new program, reduce also the value of e significantly. 

It worths mentioning that validated continuation can be applied to other delay 
equations. In particular, one interesting future project would be to apply the method 
to study periodic solutions of the Mackey- Glass equation (see [17]) 

±{t) = ^ ~ ^\.„ - Px{t), a,/3,r>0,nGN, (46) 

1 + [x[t - t)J" 

for which the existence of more than one SOPS in (46) is an open conjecture, for 
certain range of parameters. We refer to [15] for more details on this conjecture. 

The author would like to thank to Roger Nussbaum, John Mallet-Paret, Kon- 
stantin Mischaikow and Eduardo Liz for helpful discussions. Also, the author would 
like to give a special thank to Jan Bouwe van den Berg for his idea about the formu- 
lation of the bifurcation analysis presented in Section 5.4. 
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